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ABSTRACT 

Context. A highly favoured mechanism of planetesimal formation is coUisional growth. Single dust grains hit each other due to 
relative velocities caused by gas flows in the protoplanetary disc, which they follow. They stick due to van der Waals forces and form 
fluffy aggregates up to centimetre size. The mechanism of further growth is unclear since the outcome of aggregate collisions in the 
relevant velocity and size regime cannot be investigated in the laboratory under protoplanetary disc conditions. Realistic statistics of 
the result of dust aggregate collisions beyond decimetre size is missing for a deeper understanding of planetary growth. 
Aims. Joining experimental and numerical efforts we want to calibrate and validate a computer program that is capable of a correct 
simulation of the macroscopic behaviour of highly porous dust aggregates. After testing its numerical limitations thoroughly we 
will check the program especially for a realistic reproduction of the compaction, bouncing and fragmentation behaviour. This will 
demonstrate the validity of our code, which will finally be utilised to simulate dust aggregate collisions and to close the gap of 
fragmentation statistics in future work. 

Mefhods. We adopt the smooth particle hydrodynamics (SPH) numerical scheme with extensions for the simulation of solid bodies 
and a modified version of the Sirono porosity model. Experimentally measured macroscopic material properties of Si02 dust are 
implemented. By simulating three different setups we calibrate and test for the compressive strength relation (compaction experiment) 
and the bulk modulus (bouncing and fragmentation experiments). Data from experiments and simulations will be compared directly. 
Results. SPH has already proven to be a suitable tool to simulate collisions at rather high velocities. In this work we demonstrate 
that its area of application can not only be extended to low-velocity experiments and collisions. It can also be used to simulate the 
behaviour of highly porous objects in this velocity regime to a very high accuracy. A correct reproduction of density structures in 
the compaction experiment, of the coefficient of restitution in the bouncing experiment and of the fragment mass distribution in the 
fragmentation experiment show the validity and consistency of our code for the simulation of the elastic and plastic properties of 
the simulated dust aggregates. The result of this calibration process is an SPH code that can be utilised to investigate the collisional 
outcome of porous dust in the low-velocity regime. 
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1. Introduction 

In gaseous circumstellar discs, the potential birthplace of plan- 
etary systems, dust grains smaller than a micrometre grow to 
kilometre-sized planetesimals, which themselves proceed to ter- 
restrial planets and cores of giant planets by gravity-driven run- 
away accretion. Depending on their size the dust grains and 
aggregates perform motions in the disc and relative to each 
other. Brownian motion, radial drift, v ertical settling and tur- 
bulent mixing cause mutual co llisions (IWeidenschiUingl [19771 
IWeidenschilling & Cuzzilll993h . 

Since real protoplanetary dust particles are not available 
for experiments in the laboratory much of the following 
work has been carried out with dust analogues such as Si02 
(iBlum & Wurm 2008). Theoretical models also refer to micro- 
scopic and macroscopic properties of these materials. Initially 
the dust grains hit and stick on contact by van der Waals forces 



dHeim et a n [T999l) . In this proce ss, which has been investi- 
gated experimen tally dBlum et alJ 120001: iBlum & WurmI 120001: 



Krause & Blum 



ni 
1 



^ 20041) and numerically dPominik & Tielen: 

1997; Paszun & DomhiiSl2006l 12001 120091: IWada et al.ll200 ' 
2QQ8, 2009) . they form fluffy aggregates with a high degree of 
porosity. 

Due to restructuring the aggregates gain a higher mass to 
surfac e ratio and reach higher velocities. Blum & WurmI (l2000l 
l2008h and Wada et al. ( 2008) showed that collisions among them 
lead to fragmentation and mass loss. Depending on the model 
of the protoplanetary disc, which provides the kinetic colli- 
sion parameters, this means that direct gr owth ends at aggre- 
gat e sizes of a few centimetres. However, IWurm et all (I2005h 
and iTeiser & W urm (2009) have demonstrated in laboratory ex- 
periments in the centimetre regime and with low-porosity dust 
that the projectile can stick partially to the target at velocities of 
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more than 20 m s ' . Thus, colHsional growth beyond centimetre 
size seems to be possible and the exact outcome of the fragment 
distribution is crucial for the understanding of the growth mech- 
anism. 

Numerical models that try to combine elaborate pro- 
topla netary disc physics with the dust coagulation prob- 
lem (IWeid enschilling et al.' 1997; Dullemond & Dominik'2005'; 



Brauer et al. 2008; Zsom & Dullemond 2008; Ormel et al. 2007> 
2009h have to make assumptions about the outcome of collisions 



between dusty objects for all sizes and relative velocities. Since 
data for these are hardly available, in the most basic versions 
of these models perfect sticking, in more elaborate ones power- 
law fragment distributions from experiments and observations 
(iMathis et al.llT977l iDavis & Rvan|[l990t iBlum & Munci3[T99l 
are assumed. The results of alike simulations highly depends on 
the assumed fragmentation kernel. However, the given experi- 
mental references have been measured only for small aggregate 
sizes. The influence of initial parameters such as the rotation of 
the objects or their porosity have not been taken into account 
for the size regimes beyond centimetre size, although they might 
play an important role (Sirono 2004; Ormel et al. 2007)- 

A new approach to model the gro wth of protop l anetary 
dus t aggregates was re cently developed bv lGiittler et al.l(l2009ah 
and IZsom et al.l (l2009l) who directly implemented the results of 
dust aggregation experiments into a Monte-Carlo growth model. 
They found that bouncing of protoplanetary dust aggregates 
plays a major role for their evolution as it is able to inhibit fur- 
ther growth and changes their aerodynamic properties. Although 
their model relies on the most comprehensive database of dust 
aggregate properties and their collisional behaviour, they were 
unable to make direct predictions for any arbitrary set of colli- 
sion parameters and were thus obliged to perform extrapolations 
over orders of magnitude. Some of these extrapolations are based 
on physical models which need to be supported by further exper- 
iments and where this is not possible by sophisticated numerical 
models such as the one presented here. 

Sticking, bouncing, and fragmentation are the important col- 
lision outcomes which need to be implemented into a coagu- 
lation code and affect the results of these models. Thus, it is 
not only important to correctly implement the exact thresholds 
between these regimes but also details concerning the outcome 
such as the fragment size distribution, the compaction in bounc- 
ing collisions, and maybe even the shape of the aggregates after 
they merged in a sticking collision. Due to the lack of impor- 
tant input information the necessity for a systematic study of all 
relevant collision parameters arises and will be addressed in this 
work. 

Because of restrictions in size and realistic environment pa- 
rameters this task cannot be achieved in the laboratory alone. 
There has been a lot of work lately on modelling the behaviour 
of dust aggregates on the basi s of molecu lar interactio ns be- 
tween the monomers (Paszun & DominikI |2006, 2008, 120091: 
IWada et al.l l2007. 2008'. 120091) . However, simulating dust aggre- 
gates with a model based on macroscopic material properties 
such as density, porosity, bulk and shear moduli and compres- 
sive, tensile and shear strengths remains an open field since these 
quantities are rarely available. The advantage of this approach 
over the molecular dynamics method, which is computationally 
limited to a few ten thousand monomers, is the accessibility of 
aggregate siz es beyond the centi metre regime. 

Recently, lJutzi et al] (l2008h have implemented a porosity 
mo del into the smo ot h par ticle hydrodynamics (SPH) code 
by iBenz & Asphaiig ( |1994|) . It was calibrated f or pumice 
material using high-velocity impact experiments dJutzi et alJ 



I2009bh and utilised to understand the formation of an aster- 
oid family (lJutzi et aP l2009ah . However, pumice is a mate- 
rial whose strength parameters decrease when it is compacted 
(crushed). Additionally, the underlying thermodynamically en- 
hanced porosity model is designed to describe impacts of some 
kms"'. Thus, it is perfectly suitable to simulate high-velocity 
collisions of porous rock-like material. 



In contrast, collisions between pre-planetesimals occur at 
relative velocities of some tens of ms"' and compressive, 
shear and tensile streng t hs are increasing with increasing den- 
sitv dBlum & Schrapled [20041; iBlum et alj [20061; iGiittler et al.l 

Hoom 



ISchaferetan(l2007h have used an SPH code based on the 
porosity model by ISirond (l2004 to simulate collisions be- 
tween porous ice in the m s"' regime. They found that a suit- 
able choice of relations for the material parameters can pro- 
duce sticking, bouncing or fragmentation of the colliding ob- 
jects. Therefore, they stressed the importance of calibrating the 
material parameters of porous matter with laboratory measure- 
ments. Numerical molecular-dynamics simulations are about to 
use a sufficient number of monomers to be close enough to the 
continu um limit an d to provide the required material parameters 
(e.g.jPasz un & Do minik 2008, reproducing experimental results 



of iBlum & Schrapler 2004). They represent an important sup- 
port to the difficult experim ental determination of these quanti- 
ties. In lGiittler et al.l (l2009bl) we have measured the compressive 
strength relation for spherical Si02 dust aggregates for the static 
case in the laboratory and have given a prescription how to apply 
this to the dynamic case using a compaction calibration experi- 
ment and 2D simulations. We have pointed out relevant bench- 
mark features. It was shown that the code is in principal capable 
of simulating not only fragmentation but also bouncing, which 
cannot be seen in molecular-dynamics codes so far 

In this work we present our SPH code with its technical de- 
tails (Sect. |2]l, experimental reference (Sect. [3]), and numerical 
properties (Sect. H)). On the basis of the compaction calibration 
simulation we demonstrate that the results converge for increas- 
ing spatial resolution and choose a sufficient numerical resolu- 
tion (Sect. 14.21 ). We investigate the differences between 2D and 
3D numerical setup (Sect. l43b and thereby improve some draw- 
backs in Guttler et al. (2009b). Artificial viscosity will be pre- 
sented as stabilising tool for various problems in the simulations 
and its influence on the physical results will be pointed out (Sect. 
1441). 



Most prominentl y we continue the calibration process started 
in IGiittler et al.l (l200 9b) utilising two further calibration experi- 
ments for bouncing (Sect. |52]i and fragmentation (Sect. [531 ). In 
the end we possess a collection of material parameters, which is 
consistent for all benchmark experiments. Finally, the SPH code 
has gained enough reliability to be used to enhance our infor- 
mation about the underlying physics of dust aggregate collisions 
beyond the centimetre regime. In future work it will be applied 
to generate a catalogue of pre-planetesimal collisions and their 
outcome regarding all relevant parameters for planet formation. 
Jointly with experiments and coagulation models this can be im- 
plemented into protoplan etary growth simulations (iGiittler et al.l 
l2009al:IZsom etal. ] |2009l) to enhance their reliability and predic- 
tive power. 
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2. Physical Model and Numerical Method 

2. 1 . Smooth Particle Hydrodynamics 

The numerical Lagrangian particle method smoot h particle hy- 
dro dynamics (SPH) was original ly introduced by iLucvl (Il977h 
and lGingold & Monaghanl(ll977h to model compressible hydro- 
dynamic flows in a strophysical applications. La ter, the method 
has been extend ed dLiberskv & Petschekl 1990|). and improved 
extensively (e.g. Liberskv et alj|1993l: iBenz & Asphauall994t 



iRandles & Liberskvl Il996r iLiberskv et al.l Il997h to simulate 
elastic and plastic deformations of solid materials. A compre- 
hensive descriptio n of SPH and its extensions can be found in 
iMonaghanI ( 120051) . 

In the SPH scheme, continuous solid objects are discretized 
into interacting mass packages, so called "particles". These par- 
ticles form a natural frame of reference for any deformation and 
fragmentation that the solid body may undergo. All spatial field 
quantities of the object are approximated onto the particle posi- 
tions Xj by a discretized convolution with a kernel function W. 
The kernel W depends on particle distance \Xi - jc,j and has com- 
pact support, determined by the smoot hing length /; . We are 
using the standard cubic spline kernel (iMonaghan & Lattanziol 
Il985h but normalised such that its maximum extension is equal 
to one smoothing length h. 

We apply a constant smoothing length. This also allows 
to model fragmentation of solid objects in a simple way. 
Fragmentation occurs when some SPH particles within the body 
lose contact with their adjacent particles. Two fragments are 
completely separated as soon as their respective subsets of parti- 
cles have reached a distance of more than 2h so that their kernels 
do not overlap any more. 

Time evolution of the SPH particles is computed according 
to the Lagrangian form of the equations of continuum mechan- 
ics, while transferring spatial derivatives by partial integration 
onto the analytically given kernel W. 



2.2. Continuum mechanics 

A system of three partial dififerential equations forms the frame- 
work of continuum mechanics. As commonly known they follow 
from the constraints of conservation of mass, momentum and en- 
ergy. The first, accounting for the conservation of mass, is called 
continuity equation 



dp dv" 



(1) 



Following the usual notation p and v" denote density and veloc- 
ity, respectively. Greek indices run from 1 to d, the dimension 
of the problem. Einstein summing notation holds throughout the 
entire paper. In contrast to the usual SPH scheme, where the SPH 
density p, is calculated directly from the particle distribution, we 
solve the continuity equation according to 



dp/ ^-n rrij 



Pi 



dW'Kh) 



(2) 



(e.g. IRandles & Liberskvll 19961) . Here the sum runs over all in- 
teraction partners j of particle nij is the particle mass of par- 
ticle j, and denotes the kernel for the particular interaction. 
Although this approach is more expensive, as it requires to solve 
an additional ordinary differential equation for each particle, it is 
more stable for high density contrasts and it avoids artifacts due 
to smoothing at boundaries and interfaces. 



The conservation of momentum is ensured by the second 
equation 



dv'' 
"dT 



1 dcr"P 



p dxP 

In SPH formulation, the momentum equation reads 



PI Pj 



dW'Kh) 



5/ 



(3) 



(4) 



Due to the symmetry in the interaction terms, conservation of 
momentum is ensured by construction. Additionally we ap- 
ply th e standard SPH artificial viscosity (iMonaghan & Gingoldl 
Il983l) . This is essential in particular for stability at interfaces 
with highly varying densities. The influence of artificial vis- 
cosity on our simulation results is investigated thoroughly in 
Sec. 1441 

The third equation, the energy equation, is not used in our 
model. Hence, we assume that kinetic energy is mainly con- 
verted into deformation energy and energy dissipated by viscous 
effects is converted into heat and radiated away. The stress tensor 
<j can be split into a part representing the pure hydrostatic pres- 
sure p and a traceless part for the shear stresses, the so called 
deviatoric stress tensor S Hence, 



a-"" = -pd"/" + S 



(5) 



Any deformation of a solid body leads to a development of in- 
ternal stresses in a specifically material dependent manner. The 
relation between deformation and stresses is not taken into ac- 
count within the regular equations of fluid dynamics. Therefore, 
they are insufficient to describe a perfectly elastic body and have 
to be extended. The missing relations are the constitutive equa- 
tions, which depend on the strain tensor 



aP 



1 idx'" dx'!" 



2 \ dxP dx" 



(6) 



It represents the local deformation of the body. The primed co- 
ordinates denote the positions of the deformed body. 

Following Hooke's law a proportional relation between de- 
formation is assumed involving the material dependent shear 
modulus yU = i-i{p), which depends itself on the density; 



S 



ap 



(7) 



However, this is only the constitutive equation for the traceless 
shear part. For the hydrostatic part of the stress tensor we adopt 
a mo dification of th e Murnaghan equation of state which is part 
of the lSironol (l2004l) porosity model: 



Pip) = K{p'o)(p/p'o - 1), 



(8) 



where p^ is the so called reference density, the density of the 
material at zero external stress, and K{p) is the bulk modulus. 

The density dependence of the bulk K{p) and shear fi{p) 
moduU is modelled by a power law 

K(p) = 2//(p) = Koip/pd''. (9) 

Although according to ISironol (l2004h p, is the initial density of 
the material at the beginning of the simulation, we, in contrast, 
want to ensure that our dust material possesses the same bulk 
modulus K{p) even for simulations with different initial densi- 
ties. In this work the dust material has two different densities 
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at the beginning of the bouncing (Sect. I5.2ll and fragm entation 
(Sect. 1531 ) calibration setup. According to lSiron ol (12004 the ma- 
terials should feature two different p,-. As a consequence Kip), 
and in particular Kq, depends on the initial setup. Since we want 
to compare and validate Kq by using two different setups we 
have to fix a unique p, for all simulations. We choose p, such 
that K{pj) - Ko is the bulk modulus of the generic uncompressed 
dust material th at is produced by the rand om ballistic deposition 
(RBD) method ( iBlum & SchriipledllOOl . 

The time evolution of the pressure is directly given by the 
time evolution of the density (Eq. [TJ and since the pressure 
is a scalar quantity it is intrinsically invariant under rotation. 
However, in order to gain a frame invariant formulation of the 
time evolution of the deviatoric stress tensor, i.e. the stress rate, 
correction terms have to be added. A very common formulati on 
for SPH (see e.g. lBenz & Asphauglll994t lSchafer et all2007h is 
the Jaumann rate form 

rice/? / 1 \ 

2yu e"/^ - -d^^e^y 1 + 5 "^K*^ + S'^^K'", 



df 



where R"^ is the rotation rate tensor, defined by 

2 \dy^ dx" 
and e"^ denotes the strain rate tensor 



(10) 



(11) 



(12) 



To determine rotation rate and strain rate tensor and thus the evo- 
lution of the stress tensor Eq.[TO]in SPH representation, the SPH 
velocity derivatives dv"/dx^ have to be calculated. The standard 
SPH formulation, however, does not conserve angular momen- 
tum due to the discretization error by particle disorder, which 
leads to a rotational instability and in particular inhibits mod- 
elling rigid rotation of solid bodies. To avoid this error, we apply 
a correction tensor C'^ ( Schafer et al. .20071) according to 



where the correction tensor is the inverse of 
dW'j 



(13) 



til J I 

— (xf-^p- 
PJ ' 



that is 

Z.p/ J '>Z^ dx] 

This approach leads by construction to first order consistency 
where the errors due to particle disorder cancel out and the con- 
servation of angular momentum is ensured. Only this allows that 
rigid rotation can be simulated correctly. 

Finally, the sound speed of the material is given by 



(16) 



Together with Eq. |9] this relation shows that the soundspeed is 
a strong function of density. This behaviour has been seen i n 
moleculardynamics simulations by iPaszun & DominikI (l2008h . 
but there is no data available from laboratory measurements. 

Up to this point the set of equations describes a perfectly 
elastic solid body. Additionally, the material simulated in this 
work are Si02 dust aggregates, which features a high degree of 
porosity and, thus, plasticity. The modifications accounting for 
these features are described in the next section. 



2.3. Porosity and plasticity 

Following the ISironol (l2004h model, the porosity <I> is modelled 
by the density of the porous material p and the constant matrix 
density p. 



1-^ 



1 



(17) 



In the following we will use the filling factor (p - p/ps. 

We model plasticity by reducing inner stresses given by cr"^. 
For this reason we need constitutive relations describing the be- 
haviour of the material during plastic deformation. These rela- 
tions are specific for each material and have to be determined 
empirically. Particularly for highly porous materials it is ex- 
tremely diflicult to acquire them. Therefore, it is an advanta- 
geous feature of our mod e l that it is based on mea surements by 
iBlum & Schrapleii (|2004 . iBlum et al.l (l2006h and iGuttler et al.l 
(I2009bl) . 



The main idea of the adopted plasticity model is to reduce 
inner stress once the material exceeds a certain plasticity cri- 
terion. In the elastic case, described by Eqns. |7] and |8] inner 
stresses grow linearly with deformation. Hence, the material re- 
turns to its original shape at vanishing external forces. Reducing 
inner stresses, i.e. deviating from the elastic deformation path, 
reduces the internal ability of the material to restore its original 
shape. Therefore, by stress reduction deformation becomes per- 
manent , i.e. p lastic. Following and expanding the approach by 
ISironol (l2004l) . we treat plasticity for the pure hydrostatic pres- 
sure p and the deviatoric stress tensor S separately. 

For the devia t oric s tress tensor we follow th e approach by 
iBenz & Asphaud ( ll995h and lSchafer et al.1 ( l2007h assuming that 
our material is isotropic, which makes the von Mises yield cri- 
terion applicable. This criterion is characterised by the shear 
strength Y, which in our model is a composite of the compres- 
sive and tensile strengths: Y{<p) - -J'Z((t>) \T((^\. The suitabilit y 
of this choice was already demonstrated in lGiittler et al.l (l2009bh . 
Since Y{(p) is a scalar we have to derive a scalar quantity from 
S"^ for reasons of comparability. We do this by calculating its 
second irreducible invariant J2 - S'-'^S"^. Finally, the reduction 
of the deviatoric stress is implemented in the following way 



(14) s"!^ ^ fS"l^ , where / = min[y-(0)/372,l]. 



(18) 



The hydrostatic pressure is limited by the tensile strength T{^) 
for p < and by the compressive strength 2(0) for p > 0: 



P^f' \ T{cf>) 4> < 0c 



(19) 



For 4>^ < 4> ^ 4>t the material is in the elastic regime and Eq. 
[8]is applied. 0^: and (p^ denote the filling factors where the elas- 
tic path intersects the tensile strength and compressive strength, 
respectively (see Fig.lTJ. 

The pressure reduction process is implemented such that 
at any time step p is computed with Eq. [8l If for a given 
p(0) > 2(0) and 4> > <t't the pressure p(0) is reduced to 2(0). 
The deformation becomes irreversible once the new reference 
density pJ, is computed through Eq. [8] and the elastic path is 
shifted towards higher densities. Hereby also the limiting filling 
factors 0(7 and (pl are set anew. In principal there are two pos- 
sible implementations for this. (1) Plasticity becomes effective 
immediately and pJ, is computed whenever p > 2. (2) Plasticity 
becomes effective after pressure decrease, which is equivalent to 
< 0(1^. We tested both implementations. For our understanding 



R. J. Geretshauser et al.: Numerical Simulations of Highly Porous Dust Aggregates 



5 



P4< 
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Fig. 1. In the modified Sirono porosity model the regime of 
elastic deformation is limited by the compressive strength 
which represents the transition threshold to plastic compression 
for p > 2(0), and the tensile strength T{(f)), which represents the 
transition threshold to plastic tension or rupture for p < T{(p). 
Returning from one of the plastic regimes to vanishing external 
pressure via an elastic path leads to a cp'^ that differs from the 
initial (pL Hence, the material has been deformed irreversibly. 



possibility (1) is closer to the underlying physical process. In ad- 
dition it proved to be more stable. According to the benchmark 
parameters the results were equivalent. 

For the tensile regime, i.e. for (p < (p^, we do no t adopt 
the da mage and damage restoration model presented in ISirond 
(12004 . This damage model for brittle material such as rocks or 
pumi ce has been developed for SP H by Benz & Asphaug ( 19941 
Il995h and was recently used by lJutzi et alJ (l2008l l2009allb[K 
It is assumed that a material contains flaw s, which are acti - 
vated and develop u nder tensile loading ( Gr adv & Kipp|[T980l) . 
ISchafer et alJ (l2007h did not find the model applicable in their 
simulations of porous ice because it includes compressive dam- 
age effects. Brittle material like pumice and rocks tend to dis- 
integrate when they are compressed: they crush. Whereas in 
our material, highly porous SiOi dust, tensile and compressive 
strength increase with compression. This is due to the fact that 
the monomers are able to form new bonds when they get in con- 
tact. Therefore, we adopt the same approach as in the compres- 
sive regime and reduce the pressure p{(p) to T{(p) once the ten- 
sile strength is exceeded. Finally, the material can rupture due to 
plastic flow. However, material that is plastically stretched can 
be compressed again up to its full strength. By choosing this ap- 
proach a "damage restoration model" is implemented in a very 
natural way. 

Finally, a remark has to be made about energy. Apart from 
energy dissipation by numerical and artificial viscosity we as- 
sume intrinsic energy conservation. We suppose that heat pro- 
duction in the investigated physical processes is negligible. 
Therefore, our model is limited to a velocity re gime below the 
sound speed c , of the dust materi al 30 m/s, see lBlum & Wurnj 
l200i |P aszun & Dominikll2008h . By choosing the approach of 
modelling plasticity via stress reduction, we assume that most 
of the energy is dissipated by plastic deformation, since the re- 
duction of internal stresses accounts for the reduction of internal 
energy. 



Since we do not solve the energy equation thermodynami- 
cally enhanced features like any phase transition such as melting 
and freezing cannot be simulated. This scheme also does not fea- 
ture a damage model. Especially in the section about fragmenta- 
tion (Sect. 15.3b the presence of any flaws in the material cannot 
be taken into account yet, although they might have influence on 
the resulting fragment distribution. 

3. Experimental reference 

3.1. Material parameters 

The material used for the calibration exper iments are highly- 
porous dust aggregates as described by iBlum & Schraplen 
(|2004|) . consisting of spherical Si02 spheres with a diameter of 
1.5 yum. For these well defined dust aggregates it is possible 
to reproducibly measure macroscopic material parameters like 
tensile strength, compressive strength, and, potentially, also the 
shear strength as needed for the SPH porosity model (see sec- 
tion l2.3b . The tensile strength for this material was measured for 
highly porous and compa cted aggregates {cp = 0.15 .. 0.66) by 
iBlum & Schrapled (120041) . These measurements support a linear 
dependence between the tensile strength and the number of con- 
tacts per monomer (increasing with increasing <p), which yields 
the tensile strength as 



Ticp)^-{lQ'-^^'-^'^) Pa 



(20) 



The compressive strengt h was measured in t he experimental 
counterpart of this paper ( iGiittler et al.ll2009bl) with an experi- 
mental setup to determine the static omni-directional compres- 
sion (ODC), whereby the sample is enclosed from all sides and 
the pressure is constant within the sample. We found that the 
compressive strength curve can be well described by the analytic 
function 



S(^) 



(p2 ■ 



(p2-(p 



(21) 



where the free parameters were measured to be (p\ - 0.12, 
4>2 - 0.58, A = 0.58, and = 13 kPa. However, these parame- 
ters were measured with a static setup and we expect a different 
strength for a dynamic collision. The parameters (p\ and (pi deter- 
mine the range of the volume filling factor, where (p\ is defined 
by the uncompressed material and corresponds to the densest 
packing. These parameters are expected to be the same for the 
dynamic compression, while p^ and A are treated as free param- 
eters, which we will calibrate in the forthcoming sections. The 
last important material parameter to describe the dust aggregates 
is the elasticity. We can determine the Young's modulus from 
measurement and simulation of the sound speed dBlum & WurmI 
120081: iPaszun & Dominidl2008h . which is c = 30 m s ^ From 
this approach the bulk modulus for the uncompressed material 
is Kq - piC^ - 300 kPa with pj — 300 kg m ~^. However, other 
plausible calculations (iWeidling et al.ll2009l) indicate an elastic- 
ity of rather 1 kPa. We will therefore also vary this parameter. 

3.2. Calibration experiment 

As a setup for an easy a nd well-defined calibration experiment 
(see IGiittler et al.ll2009"bl) . we chose a glass bead with a diame- 
ter of 1 to 3 mm, which impacts into the dust aggregate material 
with a velocity between 0.1 and 1ms"' under vacuum condi- 
tions (pressure 0.1 mbar). We were able to measure the deceler- 
ation curve, stopping time, and intrusion depth of the glass bead 
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(for various velocities and projectile diameters) and the com- 
paction of the dust under the glass bead (for a 1 . 1 mm projectile 
with a velocity of 0.65 m s"'). These results will serve for cali- 
brating and testing the SPH code. 

For the measurement of the deceleration curve, we used an 
elongated epoxy projectile instead of the glass bead. The bottom 
shape and the mass resembled the glass bead, while the lower 
density and the therefore longer extension made it possible to 
observe the projectile during the intrusion. The projectile was 
observed by a high-speed camera (12,000 frames per second) 
and the position of the upper edge was followed with an accu- 
racy of » 3 j^m. We found that - independently of velocity and 
projectile diameter - the intrusion curve can be described by a 
sine curve 



h(t) 



-D ■ sin 



n t 
2T\ 



t < r,. 



(22) 



Here, D and Ts denote the intrusion depth and the stopping time 
of the projectile. The stopping time is defined as the time be- 
tween the first contact and the deepest intrusion at zero velocity. 
In Sect. 15.1.21 we will use the normalised form of the intrusion 
curve with D -T^ - 1 . 

For the intrusion depth we found good agreement with a lin- 
ear behaviour of 



D = 8.3 ■ 10 



-4 





m~ s 



kg 



mv 



(23) 



where v is the impact velocity of the projectile and m and 
A - nr are the mass and cross-sectional area of the projec- 
tile, respectively. We found the stopping time to be indepen- 
dent of the velocity and only depending on the projectile size 
(3.0 ±0.1 ms for 1 mm projectiles and 6.2 + 0.1 ms for 3 mm 
projectiles). 

The compaction of dust underneath the impacted glass bead 
was measured by x-ray micro-tomography. The glass bead diam- 
eter and velocity for these experiments correspond to the com- 
paction calibration setup described in Table [T] The dust sample 
with embedded glass bead was positioned onto a rotatable sam- 
ple carrier between an x-ray source and the detector During the 
rotation around 360°, 400 transmission images were taken, from 
which we computed a 3D density reconstruction with a spatial 
resolution of 21 yum. The according results of the density recon- 
struction can be found in lGiittler et alJ (l2009bl) . where we found 
that roughly one sphere volume under the glass bead is com- 
pressed to a volume filling factor of (p x 0.23, while the sur- 
rounding volume is nearly unaffected with an original volume 
filling factor of (p * 0.15. In this work, we will focus on the 
vertical density profile through the centre of the sphere and the 
compressed material (see section|4ll. 

3.3. Further benchmark experiments 

In this section, we will present two further experiments which 
will be used for the validat i on of the SPH code in sections 15.21 
and 15.31 iHeifielmann et al.l (l2007h performed low-velocity col- 
lisions (v = 0.4 m s"') between cubic-shaped, approx. 5 mm- 
sized aggregates of the material as described in Sect. 13.11 and 
found bouncing whereby approx. 95% of the energy was dis- 
sipated in a central collisi on. Detailed investiga tion of the com- 
paction in these collisions jWeidling et al.l2009l) revealed signif- 
icant compaction of the aggregates (from (f> - 0.15 to (p = 0.37) 
after approx. 1,000 collisions. The energy needed for this com- 
paction is consistent with the energy dissipation as measured by 
iHeiBelmann et alJ (l2007h . 




frogment moss / projectile moss 

Fig. 2. Cumulative mass distribution of the fragments after a dis- 
ruptive collision, which can be described by a power law. The 
divergence for low masses is due to depletion of small aggre- 
gates because of the camera resolution. 



A further experiment deals with the disruptive fragm enta- 
tion of dust aggregates (for details see lGiittler et al .112009 alibi) . In 
this case, a dust aggregate with a diameter of 0.57 mm consist- 
ing of 1.5 jum spherical SiOi dust with a volume filling factor 
of = 0.35 collides with a sohd glass target a t a velocity of 
8.4 m s" (also see Fig. 20 in lGuttler et alj2009bl) . The projectile 
fragments and the projected sizes of these fragments are mea- 
sured with a high-speed camera with a resolution of 16 fim per 
pixel. As the mass measurement is restricted to the 2D images, 
the projected area of each fragment is averaged over a sequence 
of images where it is clearly separated from others. From this 
projected area, the fragment masses are calculated with the as- 
sumptions of a spheric shape and an unchanged volume filling 
factor Fig. [2] shows the mass distribution in a cumulative plot. 
For the larger masses, which are not depleted due to the finite 
camera resolution, we find good agreement with a power-law 
distribution 



(24) 



n{m) dm — I — dm , 



where m is the normalised mass (fragment mass divided by pro- 
jectile mass) and fi - 0.22 is a measure for the strength of frag- 
mentation, being defined as the mass of the largest fragment 
divided by the mass of the projectile. The exponent a has a 
value of 0.67. A simila r distribution was already described by 
iBlum & MunciJ (Il993h for aggregate-aggregate fragmentation 
of ZrSi04 aggregates with a comparable porosity. 

4. Numerical Issues 

Before we perform the calibration process, some numerical is- 
sues have to be resolved. For instance, it is unfeasible to sim- 
ulate the dust sample, into which the glass bead is dropping in 
the compaction calibration experiment presented in section|3] as 
a whole. It is also unfeasible to carry out all necessary compu- 
tations in 3D. Therefore, we will only simulate part of the dust 
sample and explore at which size spurious boundary effects will 
emerge. Most of the calibration process has been conducted in 
2D, but the differences between 2D and 3D results will be dis- 
cussed and quantified. 
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sphere 
particles 




Fig. 3. Compaction calibration setup in 2D or, respectively, 
cross-section of 3D compaction calibration setup. A glass sphere 
impacts into a dust sample (radius rsampie) with initial velocity vq. 
The dust sample is surrounded by boundary particles at its bot- 
tom. Their acceleration is set to zero at every time step. The ver- 
tical density profile at maximum intrusion serves as calibration 
feature. 



Physical Quantity 


Symbol 


Value 


Unit 


Glass bead 


Bulk density' * 


Po 




kgm ' 


Bulk modulus'** 


Ko 


5 X 10' 


Pa 


Mumaghan exponent'** 


n 


4 


- 


Radius 


r 


0.55 X 10"' 


m 


Impact velocity 


vo 


0.65 


ms"' 


Dust sample 








Initial density 


Pi 


300 


— 1 —I — 

kgm 


Bulk density 


P.< 


2000 


kg m"' 


Reference density 


Po 


300 


kgm~' 


Initial filling factor 


<p, 


0.15 




Bulk modulus 


Ko 


3x 10' 


Pa 


ODC mean pressure 


Pm 


1300 


Pa 


ODC max. filling factor 


<Pl 


0.58 




ODC min. filling factor 


01 


0.12 




ODC slope 


A 


0.58 





Table 1. Numerical parameters for the compaction calibration 
setup. ODC stands for omni-directional compression relation 
(Eq.lTTTi. Qu antities marked by (*) represent the parameters for 
sandstone in iMeloshlJTosl which we adopt for glass here. 



In this context we make use of 2D simulations in cartesian 
coordinates although the symmetry of the problem suggests the 
usage of cylindrical coordinates. However, the SPH scheme in 
cylindrical or polar coordinates battles with the problem of a sin- 
gularity at the origin of the kernel fu nction. There e xist only few 
attempts to solve this issue (e.g. Om ang et al.ll2006l) . but they are 
still under development and require high implementation efforts. 
Since in our case 2D simulations only serve as indicator for cal- 
ibration and 3D simulations are aimed at, we stick to cartesian 
coordinates. 

The glass bead is simulated with the Murnaghan equation of 
state 



P(P)= - 
\ n 



1 



(25) 



following the usual laws of continuum mechanics as presented 
in section 12.21 The compaction calibration setup is initialised 
with the numerical parameters shown in Table [T] unless stated 
otherwise in the text. Our tests showed that the maximum in- 
trusion depth and the density profile are the calibration parame- 
ters which are most sensitive to changes of the numerical setup. 
Density profiles (e.g. Figs. |4] and |5j display the filling factor 
along a line through the centre of the sphere and perpendicular 
to the bottom of the dust sample (see Fig. |3]l. The origin in the 
diagrams represents the surface of the unprocessed dust sample. 
The glass sphere has been removed in these figures. 

4.1. Computational domain and boundary conditions 

In 2D simulations we tested the effect of changing size and shape 
of the dust sample. Initially the particles were set on a triangular 
lattice with a lattice constant of 25 /im. To be geometrically con- 
sistent with the cylindric experimental setup, firstly, we utilised 
a box (width 8 mm) and varied its depth: 1.375 mm, 2.2 mm, 
3.3 mm, and 5.5 mm. This is equivalent to 2.5 x,4x,6x, and 
10 X Tsphere- Comparing the density profiles (Fig. |4l top), two 
features are remarkable: (1) The maximum filling factor at the 
top of the dust sample {<p ^ 0.27 at D ^ -0.6 mm) and the in- 
trusion depth D is nearly the same for all dust sample sizes. (2) 



For lisampie < 3.3 mm we find spurious density peaks at the lower 
boundaries (D » -1.4mm and D x -2.2 mm). 

In order to reduce the computation time we simulated the 
dust sample as a semicircle with the same radius variation as 
above. The resulting density profiles are shown in Fig. |4] (bot- 
tom). In contrast to the corresponding simulations with the box- 
shaped samples we find for r^ampie ^ 1 -375 mm an increased 
maximum filling factor and a slightly reduced intrusion depth. 
Due to the greater amount of volume lateral to the intrusion 
channel, material can be pushed aside more easily than inside the 
narrow boundaries of the semicircle. Therefore, a higher frac- 
tion of the material is compressed to higher filling factors. For 
'"sample > 3.3 mm the spurious boundary effects become negli- 
gible within the compaction calibration setup and the density 
structure shows no significant difference for box-shaped and 
semicircle shaped dust samples. 

Hence, all computations of section |4] are conducted on the 
basis of a semicircle in 2D or a hemisphere in 3D with a radius 
of ^sphere = 3.3 mm. 

In all cases the dust sample is bordered by a few layers of 
boundary particles. The acceleration of these particles is set to 
zero at each integration step, simulating reflecting boundary con- 
ditions. Apart from that, i.e. in terms of equation of state, they are 
treated like dust particles. We also tested damping boundary con- 
ditions by simulating two layers of boundaries. The outer layer 
was treated as described above, the inner (sufficiently large) 
layer was simulated with a high artificial a-viscosity. Since there 
was no significant difference in the outcome we fix all bound- 
aries in the afore mentioned way and treat them as reflecting. 

4.2. Resolution and Convergence 

LTutzil (I2008h found in his studies of a basalt sphere impacting 
into a porous target that the outcome of his simulations strongly 
depends on the resoluti on. With a calibration setup similar to the 
one used in this paper iGeretshauseJ (l2006h confirms that also 
in simulations of the type presented in this work a strong res- 
olution dependence is present. He has found that the intrusion 
depth of the glass bead can be doubled by doubling the resolu- 
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-4 -3 -2 -1 

sample height [mm] 

Fig. 4. Vertical density profile at maximum intrusion for the 
compaction calibration setup and diff'erent shapes of the 2D dust 
sample (box and semicircle). Depth dsampie of an 8 mm wide box 
(top) and radius rsampie of the semicircle (bottom) were varied. In 
both cases spurious boundary effects appear for t/sampie < 3.3 mm 
and Tsampie < 3.3 mm, respectively. 



tion. Sin ce the calibration experiments presented in lGiittler et al.l 
(l2009bl) are extremely sensitive even to minor changes of the 
setup, the convergence properties of porosity model and underly- 
ing SPH method will be investigated carefully in this paragraph. 
Additionally, we will study the differences in the outcome of 2D 
and 3D setup. 

For the 2D convergence study particles were initially put on a 
triangular lattice again. The lattice constants Ic were 100, 50, 25, 
and 12.5 fim for the compaction calibration setup. The smooth- 
ing length h was kept constant relative to at a ratio of 5.6 x Ic- 
The maximum number of interaction partners was /^ax ~ 180, 
the average /av ~ 100 and the minimum I^in ~ 30. 

In the 3D convergence study we are using a cubic lattice with 
edge lengths 1^ = 100, 50, and 25 yum. The latter was simulated 
with 3.7 million SPH particles, which represent the limit of our 
computational resources. We fixed h - 3.75 x Ic which yielded 
^max ~ 370, /av ~ 240, and ~ 70. The results are presented 
in Fig.|5] In contrast to the plots in Fig.|4]the glass sphere is not 
removed here. Coming from the right side of the plot, the filling 
factor rapidly decreases from a high value outside the scope of 
the plot indicating the sphere. The filling factor reaches its min- 
imum at an artificial gap between sphere and surface of the dust 
sample. The width of this gap is about one smoothing length 
h. The existence of the gap has two reasons: (1) Sphere mate- 



-e- 
o 

s 



-e- 

u 
O 
tj 




-2 -1.5 -1 -0.5 

sample height [mm] 

Fig. 5. Convergence study of the vertical density profile for 2D 
(top) and 3D (bottom) compaction calibration setups. The filling 
factor increase towards the surface of the dust sample accounts 
for the glass bead which is not removed in this plot. Simulations 
were performed for different spatial resolutions. All curves show 
a characteristic filling factor minimum between sphere and dust 
sample and a characteristic filhng factor maximum indicating 
the dust sample surface. 



rial and dust material have to be separated by artificial viscosity 
due to stability reasons. This will be discussed below. (2) The 
volume of the sphere represents an area of extremely high den- 
sity and pressure with respect to the dust sample. This area is 
smoothed out due to the SPH method. The width of the smooth- 
ing is given by the smoothing length. Although a clear conver- 
gence behaviour is visible in Fig. |5] for both the 2D and the 
3D case, a more unique convergence criterion has to be found. 
For this purpose we choose the maximum intrusion depth which 
proved to be very sensitive to resolution changes. The shape of 
the filling factor profile offers two choices to determine the intru- 
sion depth: (1) The filling factor minimum which represents the 
middle between sphere and dust sample and (2) the filling fac- 
tor maximum (peak) of the dust material left to the gap between 
sphere and dust sample. 

Fig. |6] shows the results for both cases in 2D (top) and 3D 
(bottom). The error bars around the minimum values represent 
the smoothing length and give an indication of the maximum er- 
ror The position of the density peak remains almost constant, 
converging to a; -0.9 mm (2D) and x -0.65 mm (3D), respec- 
tively, for higher resolutions. The position of the density min- 
imum at low resolutions significantly differs from the position 
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Fig. 6. Convergence study of the maximum intrusion depth for 
the 2D (top) and 3D (bottom) compaction calibration setups. 
Filled symbols represent the position of the filling factor peak 
of the dust material, whereas empty symbols denote the position 
of the minimum filling factor at the gap between glass bead and 
dust material. The values are derived from the density profiles 
in Fig. |5] The smoothing length is indicated by the error bars. 
While the peak position remains almost constant at ^ -0.9 mm 
(2D) and ^ -0.65 mm (3D) with increasing spatial resolution, 
the position of the filling factor minimum quickly converges to 
the same value. This is due to the artificial separation of dust and 
glass materials, which is in the order of a smoothing length. 



of the density peak, but converges quickly to the same intrusion 
depth with higher resolutions. However, the differences between 
the extrema remain well within one smoothing length. This is 
due to the fact that sphere and dust sample are separated by about 
one smoothing length. Comparing 2D and 3D convergence the 
3D case seems to converge more quickly. 

Due to the findings of this study we choose a spatial resolu- 
tion of Ic = 25 fim for further simulations in 2D. In the 3D case 
Ic - 50 jum is sufficient, but Ic < 50 /urn is desirable if it is feasi- 
ble. After defining suitable values for the spatial resolution we 
now turn to the numerical resolution, which for the SPH scheme 
is given by the number of interaction partners of each single par- 
ticle. For the investigation of this feature we performed a study 
utilising the 2D compaction calibration setup with a spatial res- 
olution of Ic - 25 fim and varied the ratio between smoothing 
length and lattice constant h/lc from 2 to 7 in steps of one. 
h/lc determines the initial number of interaction partners that is 
smoothed over. The resulting maximum, average, and minimum 
interactions /max, hy, and and the corresponding smoothing 
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Fig. 7. Convergence study for the density profile using the 2D 
setup and varying the smoothing length /?. Through this variation 
the number of interaction partners is varied according to Table|2] 
The glass bead has been removed in this plot. For h < 0.075 mm 
clear signs of instabilities are visible. For /i > 0.1 mm the filling 
factor has the same value and its position remains constant. The 
smoothing of the boundary of the dust sample is increased for 
increasing h. 

Table 2. Parameters for the convergence study regarding inter- 
action numbers 



h 


h/lc 




/av 




T 




0.050 mm 


2 


3 


13 


25 


16.2 h 


132401 


0.075 mm 


3 


10 


30 


55 


14.8 h 


96912 


0.100 mm 


4 


16 


53 


92 


19.8 h 


71175 


0.125 mm 


5 


23 


82 


142 


19.0 h 


56347 


0.150 mm 


6 


32 


116 


205 


21.6 h 


46782 


0.175 mm 


7 


43 


158 


274 


24.3 h 


39980 



lengths h can be found in Table|2] Additionally, we measured the 
computation time T^omp that the simulations took on 4 cores of a 
cluster with Intel Xenon Quad-Core processors (2.66 GHz) for a 
simulated time of 5 ms and the number of integration steps A^step 
of our adaptive Runge-Kutta Cash-Karp integrator 

Comparing the density profiles in Fig. Q (where the glass 
bead has been removed) instabilities in the form of filling fac- 
tor fluctuations due to insuflicient interaction numbers appear 
for smoothing lengths h < 0.075 mm, i.e. for /av < 30. For 
h > 0.1 mm the density profile maintains essentially the same 
shape; The position and height of the filling factor peak remains 
nearly the same and ^ smoothly drops to ^ 0.18 towards the bot- 
tom of the dust sample. Only the sharp edge at the top of the dust 
sample is smoothed out over a wider range due to the increased 
smoothing length. 

Table|2]shows that the number of integration steps A^step is de- 
creasing with increasing interaction numbers. This is due to the 
fact that elastic waves inside the dust sample are smoothed out 
over a wider range causing the adaptive integrator to increase the 
duration of a time step, since density fluctuations do not have to 
be resolved as sharply as at smaller smoothing. As expected, the 
computation time Tcomp is generally increasing with increasing 
number of interactions. There are two exceptions; h - 0.075 mm 
and h - 0.125 mm. Here, the decrease of A^steps overcompen- 
sates the increase of interactions leading to a decrease of T^omp- 
Hence, a ratio h/lc ~ 5 yields the necessary accuracy and an ac- 
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ceptable amount of comput ation time. Th i s study also justifies 
the choice of h/lc - 5.6 in iGiittler et al.l (l2009bh and we will 
stick to this ratio throughout this paper. 

According to these findings, for 3D simulations theoretically 
an average interaction number of lH"^ x 750 would be needed 
to achieve the same numerical resolution. However, alike sim- 
ulations are unfeasible and our choice of I^y x 240 in 3D is 
equivalent to I^y ^ 40 in 2D which should provide sufficient and 
rehable accuracy. 



4.3. Geometrical difference - 2D and 3D setups 

As one can easily see in Fig.|6]2D and 3D simulations have sig- 
nificantly different convergence values for the intrusion depth. 
This deviation is caused by the geometrical difference of the 2D 
and 3D setup. The 2D setup (glass circle impacts into dust semi- 
circle) represents a slice through a glass cylinder and a semi- 
cylindrical dust sample, which implies an infinite expansion into 
the third spatial direction. Whereas the 3D setup represents a real 
sphere dropping int o a "bowl" of dust. Th e relation for the intru- 
sion depth found bv IGiittler et al.l (l2009bl) (see Sect. [3]) contains 
a geometrical dependence: D oc mv/A, where D is the intrusion 
depth, m the mass of the impacting glass bead, v its impact ve- 
locity, and A its cro ss-section. 11120 is a mass per unit length. 
IGiittler et al.l (l2009t)l) already exploited this relation in order to 
determine a rough correction factor between 2D and 3D simual- 
tional setups: 



is nr~p ■ V 
3^ 2r 



37T A2D 



(26) 



Hence, the 2D intrusion depth has to be corrected by a factor of 
» ^ in order to determine the 3D intrusion depth. The compar- 
ison is shown in Fig. |8] Again we choose the 2D and 3D data 
gained in the convergence study for the peak filling factor values 
shown in Fig.|6] Fig.|8]shows the original 2D data, the corrected 
2D data, and the according 3D data (with error bars displaying 
the smoothing length). The 3D values nicely follow the rough 
correction and remain well within the maximum er ror. This 
compar ison justifies the correction of the results in Giit tler et al.l 
(l2009t)l) . With the aid of this - now verified - correction factor 
the 2D data of the intrusion depth can be converted into 3D data 
and all calibration tests involving the intrusion depth can be car- 
ried out in 2D, which saves a significant amount of computation 
time. Comparing the vertical density profiles of 2D and 3D se- 
tu ns in Fig.|5]resolves a lso another issue that remained unsolved 
in IGiittler et al.l ( l2009t)l) . According to the experimental data the 
filling factor drops to a value of <p x 0.16 within k 0.6 mm away 
from the bottom point of the glass bead. For high-resolution 2D 
simulations (Fig. jSj top) the filling factor does not drop to this 
value within the whole dust sample. However, the 3D simula- 
tions (Fig. |5] bottom) show that this effect again is due to the 
difference of 2D and 3D geometry. With the 3D setup the fill- 
ing factor drops to a; 0.16 within a; 0.9 mm. All deviations 
from experimental findings due to this effect, in particular the 
presence of a large amount of volume with (p < 0.2 in the cumu- 
lated v olume over filling factor diagram (Fig. 15 in IGiittler et al.l 
I2009bl) . can in principal be removed by switching to 3D simula- 
tions. They do not represent a fundamental error in the porosity 
model. 
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Fig. 8. Verification of the 2D-3D correction factor. Filled sym- 
bols denote the position of the filling factor peak of the dust 
material in Fig. |5] Triangles represent 3D and squares 2D val- 
ues. The conversion from 2D to 3D intrusion depth utilising the 
correction factor from Eq. |26]and Eq. |23]due to the geometrical 
difference is indicated by the line without symbols. The 3D val- 
ues are in very good agreement with the very rough theoretical 
prediction. They lie well within the errors. 



4.4. Artificial Viscosity 

Since artificial viscosity plays an eminent role for the stability 
of SPH simulations we investigate its influence on the outcome 
of our compaction calibration setup. Only artificial (j-viscosity 
is applied in our test case, but in a threefold way: (1) It damps 
high oscillation modes of the glass bead caused by the stiff 
Murnaghan equation of state (EOS). Thereby it enlarges the time 
step of our adaptive integrator and saves computation time. (2) 
It is used to provide the dust material with a basic stability. (3) 
It separates the areas of Murnaghan EOS and dust EOS and pre- 
vents a so called "cannonball instability". For all three cases also 
the influence of yS-artificial viscosity has been tested, but its in- 
fluence on all benchmark parameters was negligible. 

(1) The choice of the a-viscosity of the glass bead proved 
to be very uncritical. We choose the canonical value a - 1.0. 
There was no influence on the physical benchmark parameters 
for all a values, except a = which produces an instability. 
Values for a > 1.0 show no significant effect on the damping 
and the influence of 0.1 < a < 1.0 on it is not too high, but still 
observ able. Hence, w e stick to the canonical value. 

(2) ISironol(l2004l) applies no artificial viscosity to his porous 
ice material because of its spurious dissipative properties. Our 
findings, shown in Fig. |9l confirm this assessment regarding 
the dust material. Within our 2D compaction calibration setup 
{Ic = 25 jum, h/lc = 5.6) we vary a from to 2 and observe its 
influence on the density profile (Fig. |9] top) and the maximum 
intrusion represented by the filling factor peak of the dust mate- 
rial (Fig.|9] bottom). 

The position of the filling factor peak ranges from a 
-0.92 mm with a = 0.0 to a; -0.62 mm at or = 2.0. This clearly 
demonstrates the dissipative feature of the a-viscosity, since a 
smaller amount of kinetic energy of the glass bead is transformed 
into plastic deformation with higher a. The residual energy must 
have been dissipated. However, the cf-viscosity-intrusion curve 
seems to saturate at a value of a -0.6 mm. The decrease of 
the maximum intrusion can also be seen in the density profile 
(Fig- 12 top). While the profile maintains nearly the same shape. 
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the height of the filling factor peak is decreasing with increas- 
ing a. Hence, an increasing artificial viscosity diminishes the 
peak pressure during compaction, which is via the compressive 
strength relation E(0) directly responsible for the height of the 
filling factor peak^ 

In contrast to ISirond (l2004 we find it necessary to apply a 
small amount of cc- viscosity to the dust material. For ff < 0. 1 the 
results show traces of an instability, which is also responsible for 
a rapid increase of the maximum intrusion. Therefore, we find it 
convenient to apply an artificial viscosity with a = 0. 1 to the dust 
material, which holds for the previous simulations of this section 
as well as the following. The choice of a non-zero a, however, 
is also justified by experimental findings: After impacting into 
the dust sample the glass bead shortly oscillates due to the elas- 
tic properties of the dust. This oscillation is damped by internal 
friction, which we model with artificial viscosity. Therefore, by 
choosing a non-zero a we take into account the dissipative prop- 
erties that our dust material naturally has. A quantitative calibra- 
tion of this parameter, however, has to be left to future work. 

(3) During our first simulations with the 2D compaction cali- 
bration setup we observed what is sometimes described in the lit- 
erature as "cannonball instability": During the compaction pro- 
cess, when the glass bead intrudes into the dust material, single 
particles at the sphere's surface start to oscillate between the do- 
mains of the Murnaghan EOS and dust EOS. Due to the severe 
discrepancy of the "stiffness" of these two equations of state the 
particles collect a huge amount of kinetic energy until they are 
fast enough to generate a pressure on the dust material that ex- 
ceeds the compressive strength 2(0). Eventually they disengage 
from the sphere's surface like a cannonball and dig themselves 
into the dust sample causing a huge amount of unphysical com- 
paction. We tackle this problem by applying to all SPH particles 
with dust EOS, which interact with glass bead SPH particles, 
the same amount of a-viscosity as to the sphere, i.e. a - 1.0. 
In our simulations this is sufficient to prevent the "cannonball 
instability". The spurious dissipation caused by this measure is 
negligible. 
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Fig. 9. Density profile (top) and maximum intrusion (bottom) for 
different values of artificial ff-viscosity. The shape of the density 
profile hardly changes, but increasing a-viscosity decreases the 
maximum filling factor and the maximum intrusion depth. 



5. Calibration 

5. ) . Compressive Strength - Compaction Properties 

In this section we will refine and extend a study of the com- 
paction properties of the dust sample, whi ch was already carrie d 
out in a similar, but less detailed way in iGiittler et al.l (l2009bl) . 
There it turned out that the quantity having the most influence 
on the compaction was the compressive strength relation Z(0) 
(Eq. I2TI ). which was measured in an omni-directional and static 
manner In order to adopt this relation for dynamic compression 
the mean pressure pm and the "slope" of the Fermi-shaped curve 
A can be treated as free parameters. The upper and lower bound- 
aries of the filling factor <p2 and 0i, respe c tively, r emain constant 
even in the dynamic case. IGiittler et al.l (l2009bl) found that by 
lowering most of the features of the compaction calibration 
setup can be reproduced in a very satisfactory manner. Hereby 
the filling factor over compressive strength curve (see also Fig. 
2 in Guttler et al. 2009b) is shifted towards lower pressures and 
the yield pressure for compression is lowered. Using only 2D 
simulations and a rough parameter grid lGiittler et al.l ( 2009bh fix 
p„, = 1.3 kPa. The "slope" A has not been considered. In this 
work we will consider A and we will perform more accurate pa- 
rameter studies for p,„. From the latter we will predict a rea- 
sonable choice for p,,,, which will represent the basis for a 3D 
simulation of the compaction calibration setup. The results of 



this simulation will be compared to results from t he laboratory. 
For the comparison we use the same features as IGiittler et al.l 
(l2009bl) . 



5.1 .1 . Fixing free parameters 

Since there are no empirical data available for p,„ in the dynami- 
cal compressive strength curve we perform a parameter study in 
order to determine a suitable choice for this important quantity. 
For this study we make use of the 2D compaction calibration 
setup and vary p„, from 0.13 to 13.0 kPa (Fig.[TQland[TTI). where 
13.0kPa represents the value for the static compressive strength 
curve. The effect of lowering p,„ can most clearly be seen in 
the vertical filling factor profile (Fig.[TOl top). First of all, more 
material can be compressed and it is compressed to higher filling 
factors. As a consequence the glass bead intrudes deeper into the 
dust sample. From experimental results we expect an intrusion 
depth of about one sphere diameter (» 1 mm). With the aid of 
the empirical relation between momentum over impactor cross- 
section ;7zvA ' and intrusion depth D (Eq.l23Tl as well as the cor- 
rection factor between 2D and 3D intrusion depth (Eq. l26b we 
estimate that D^^ « 1 mm corresponds to D'^ x 1.42 mm. Fig. 
[TT]shows the maximum intrusion over the stopping time for var- 
ious values of p,,, (labels). The estimated D^^ is indicated by a 
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Fig. 10. Density profile (top) and maximum peak filling factor 
from the density profile over p,„ (bottom) for different values 
of pi„, which represents the mean pressure in the compressive 
strength relation 2(0) (Eq. 1211 1. Lowering p,„ from 13.0 kPa 
(static compressive strength) increases the compressed volume, 
its filling factor, and the maximum intrusion depth of the glass 
bead. The parameter study has been performed using the 2D 
compaction calibration setup. 



dashed line. We deduce that regarding intrusion depth a dynamic 
mean pressure pm - 0.26 kPa is a suitable choice. 

This is supported by the peak filling factor appearing in the 
filling factor profile (Fig. (TO] bottom). Empirical data indicate 
that in case of the compaction calibration setup a peak filling 
factor 0n,ax ~ 0.3 can be expected. The comparison between 2D 
and 3D results (Sect. 14.31 ) has shown that the peak filling factor 
in the vertical density profile in the 2D case is generally higher 
than for the same situation in 3D. The equivalent of (p^^^ x 0.3 
is a maximum filling factor of (p~^^^ ^ 0.34 in 2D. This points to 
a choice of p,„ ^ 0.3 kPa, which is consistent with the findings 
from the intrusion depth. 

The compressive strength relation S((^) (Eq. l2Tl i contains a 
second free parameter A, which accounts fo r the "slope" of the 
Fermi-shaped curve. In lGiittler et alJ (l2009bh this parameter was 
chosen to be the same as the one of the static omni-directional 
compressive strength curve and a closer investigation was not 
carried out. In order to understand its effect on the compaction 
properties of the dust sample we are utilising the 2D compaction 
calibration setup again and vary A from 0.55 to 0.80. The results 
are presented in Fig. [12] From the vertical density profile (Fig. 
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Fig. 11. Maximum intrusion over stopping time for different val- 
ues of the mean pressure p,„ (labels, in kPa) in the compressive 
strength relation 2(0) (Eq.lSTb using the 2D compaction calibra- 
tion setup. The dashed line indicates the 2D intrusion depth that 
is equivalent to a 3D intrusion depth of 1 mm according to 
Eq. |23]and the 2D-3D correction factor from Eq. [26] This sup- 
ports the choice of the mean pressure p„, = 0.26 kPa for further 
simulations. 



[T2I top) it can be seen that increasing A increases the intrusion 
depth, but not as effectively as by lowering the mean pressure 
p,„. This is due to the fact that increasing A is hardly increasing 
the peak filling factor in the vertical density profile. More vol- 
ume is compacted to lower filling factors. Whereas by lowering 
p„i the total amount of volume that is compacted is smaller, but 
it is compacted to higher filling factors. This behaviour can be 
seen comparing the density profiles for the p„, variation (Fig.[TOl 
top) and A variation (Fig.[T2| top) and particularly from the cu- 
mulated volume over filling factor curve (Fig. [12] bottom). This 
figure shows the normalised fraction of compacted volume cor- 
responding to a volume filling factor > (p. The intersection of the 
curves with the y-axis represents the total compressed volume, 
which is increased from x 7 to x 9.5 sphere volumes. According 
experimental measurements here expect a value of roughly one 
sphere volume (see section l3?2l i. Especially for 0.18 < (p < 0.23 
the compacted volume fraction is increase d. The equivalent fig- 
ure for the p,„ variation can be found in iGtittler et al.l (l2009bl 
Fig. 15). From the comparison of 2D and 3D calibration setups 
(section 14.3b we know that a huge amount of this compaction 
(especially in the lower filling factor regime) is due to the ge- 
ometrical difference. However, the experimental data do not in- 
dicate a particularly high amount of compaction to lower filling 
factors (rather the contrary) and, therefore, we maintain our ini- 
tial choice of A = 0.58. 



5.1.2. Comparison with experiments 

For the comparison with experiments we performed a simulation 
using the 3D compaction calibration setup (see Table [TJ with 
two exceptions: (1) The bulk modulus of the dust material was 
set to Kq - 2 kPa (instead of Kq - 300 kPa) since findings pre- 
sented below indicated a much lower bulk modulus. However, 
the choice of Ko has little influence on the compaction proper- 
ties calibrated for in the compaction calibration setup. It is more 
important for bouncing and fragmentation, which will be shown 
below. (2) The mean pressure of the compressive strength re- 
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Fig. 12. Density profile (top) and normalised volume fraction of 
compacted volume over its corresponding filling factor (bottom) 
for different values of A in the compressive strength relation (Eq. 
I2TI ). Increasing A increases the amount of compacted volume 
at filling factors 0.18 < (f> < 0.23 and thereby the maximum 
intrusion depth. Whereas the peak filling factor in the density 
profile is not changed. 



lation (Eq. |2TJ was set to p,„ = 0.26 kPa as suggested in the 
previous section. 

The following features of the compaction calibration setup 
were measured in the laboratory and will be used here for com- 
parison: (1) the stopping time Ts, (2) the intrusion curve of the 
projectile, (3) the vertical density profile, (4) the cumulated vol- 
ume over filling factor curve, and (5) the cross-section through 
sphere and dust sample displaying the filling factor 

(1) The experiments show that the stopping time T, of the 
glass bead, i.e. the time from first contact with the dust sample 
until deepest intrusion, is nearly constant at r,"'' - 3.0 + 0.1 ms 
for 1 mm projectiles over different impact velocities (see section 
13. 11 1. In our simulation we find Tf™ - 2.42 + 0.05 ms, which is 
not in excellent agreement but also not too far off the experimen- 
tal results. 

(2) The intrusion curve /i(f) was cleared from gravity effects 
and normalised through h'{t') - h{t)/D and f' - t/Ts, where 
h{t) is the position of the bottom of the glass bead as a function 
of time, D the maximum intrusion depth, and the stopping 
time. At first contact t is h'jf = 0) = and at deepest intrusion 
h'(t' = 1) = -1 (see also lGiittler et al.l2009bi section 3.2.2). The 
comparison is shown in Fig. [T3] The intrusion curve generated 



by our simulation lies well within the data from the experiments 
with 1 mm and 3 mm spheres and only slightly below the sine 
curve fitted to the empirical data. 

(3) The comparison between the experimentally measured 
vertical density profile and the result from our simulation is 
shown in Fig. [14] The crosses represent the data from two exper- 
iments where the sphere has not been removed for the measure- 
ment. For this reason the filling factor reaches extremely high 
values at X! - 1 mm indicating the bottom of the glass bead. The 
vertical density profile of our simulation is given by the solid 
line and the vertical dashed line is placed at its filling factor 
peak at -1.02 mm representing the maximum intrusion depth. 
Comparing with the experimentally measured maximum intru- 
sion depth of -1 .07 mm this is an excellent result. Since a depth 
of X! -1 mm was aimed for using Eqns. [23] and [26l and the 2D 
intrusion depth study of the previous section this result also sup- 
ports the validity of these relations. In addition to the exact value 
of the intrusion depth our simulation also reproduces the shape 
of the given experimental vertical density profile very well. Only 
the step-like structure at -1.5 mm does not find an equivalent 
in our simulation but it is nicely interpolated. 

(4) The vertical density profile shows only a cut through 
the compressed volume. It contains information about the exact 
structure of the compression. The cumulated volume over fill- 
ing factor curve (Fig.[T5]l has the advantage that it represents the 
total compressed volume together with its filling factors, hence, 
these features are not fully independent but focus on different as- 
pects of the compression. The cumulated volume is normalised 
through the sphere volume. It accounts for the volume fraction 
of compacted area corresponding to a volume filling factor > 0. 
In general, the experimental reference and our simulation show 
a very good agreement. Slightly too much volume is compressed 
to high filling factors in the simulation which leads to an almost 
constant deviation for <p < 0.26. However, the slope is repro- 
duced very well. 

(5) The comparison between the cross-sections through 
sphere and dust sample along the z-axis (Fig.[T6]l reveals reasons 
for the excess of compressed volume. First, the cross-section of 
the sphere is artificially enhanced by the smoothing of its bound- 
aries, which is inherent of the SPH method. One effect of the 
smoothing is the existence of a gap between sphere and dust 
sample, which was already discussed in Sect. 14.21 and is also 
clearly visible in Fig. [16] (top). Hence, we assume that the dust 
sample actually begins, where it has its maximum compression. 
The fact that the sphere pokes out of the dust sample a bit more 
than in the experiment has its reason also in the artificial en- 
largement of the cross-section. Second, it can be seen that in the 
experimental reference (Fig. [16] bottom) the compacted region 
is much narrower and more concentrated underneath the sphere. 
In the simulated result the compacted region is a bit broader 
This indicates that the shear strength seems to be a bit smaller 
than we assume. Third, the compaction reaches too high filling 
factors, which was already visible in the cumulated volume di- 
agram. Nevertheless, both cross-sections match very well, es- 
pecially with respect to the mediocre resolution. Remarkably, 
even the slight intrusion channel on the left and right side of the 
sphere, which features a slight compression, can be reproduced. 

5.2. Bulk Modulus 

As it was mentioned in section 13.11 there are two estimates 
for the bulk modulus Kq for the uncompressed material with 
<p w 0.15. A^o in our model is also the pre-factor for comput- 
ing the bulk modulus for higher filling factors with the aid of 
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(a) (b) (c) (d) 

Fig. 17. Bouncing sequence for f = ms (a), t - 10 ms (b), t - 18 ms (c), and f = 25 ms (d). The colour code indicates the 
filling factor. An aggregate consisting of dust particles (Sirono EOS, diameter 1.0 mm) hits a solid surface simulated by boundary 
glass particles (Murnaghan EOS, diameter 1.6 mm, thickness 0.1 mm) with a velocity of vq = 0.2 ms"'. For this simulation a bulk 
modulus of Kq =5.0 kPa and p„i - 0.26 kPa were used. The aggregate hits the surface and starts to be compacted at its bottom 
(b). While the plastic deformation at the bottom increases, the aggregate is also deformed elastically: it gets broader (c). Eventually 
it leaves the surface with a final velocity v/ (d). It features a permanent compaction, while the elastic deformation vanishes. (An 
animation of this figure is available in the online journal.) 
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Fig. 13. Comparison between intrusion curves from drop exper- 
iments using 1 mm and 3 mm spheres and our 3D simulation 
(mean pressure p,„ = 0.26 kPa). The curves are normalised such 
that the origin represents first touch of sphere and dust sam- 
ple and (1,-1) denotes stopping time at maximum intrusion. The 
simulated curve lies slightly underneath the fit to the experimen- 
tal data, but well within the errors. The deviation is due to a 
smaller stopping time than in the experiments. 





0.45 




0.4 




0.35 


-e- 




u 


0.3 


O 




tj 


0.25 








0.2 








0.15 








0.1 




0.05 










Experiment 1 * i 




■ Experiment 2 *i 




p,^ - 0.26 kPa i 




«i 





















-3 -2.5 -2 -1.5 -1 -0.5 

sample height [mm] 

Fig. 14. Comparison between experimentally measured (crosses, 
sphere not removed) and simulated (solid line, sphere removed) 
density profiles at maximum intrusion for the compaction cali- 
bration setup. The dashed line indicates the position of the sim- 
ulated maximum intrusion depth given by the density peak at 
-1.02 mm. The simulation was carried out in 3D using a mean 
pressure = 0.26 kPa for the compressive strength relation 
(Eq.l2ni. Both profiles are in excellent agreement. The fact that 
the step-like structure of the experimental data cannot be seen in 
the simulation is a minor drawback since it is interpolated nicely. 



a power law (Eg. [9| l. Since these indirect findings, 1 kPa (see 
IWeidling et al. 2009") and 3 00 kPa (from sound sp eed measured 
bv lBlum & Wurmi 2008: Pa szun & Dominikl2008h . differ by two 
orders of magnitude, we try to find a suitable value for utilis- 
ing a calibration experiment. 

Sim ulating the low velocity collision setup bv lWeidling et alJ 
(l2009l) . a 3D dust sphere (0, = 0.15,1 mm diameter, 267737 par- 
ticles, Ic - 12.5 yum) drops onto a solid surface (cylindrical, 
1.6 mm diameter, 0.1 mm thick, 1 15677 particles, = 12.5/Ym) 
with initial velocity vq = 0.2 ms"' (see Fig. [TtTi. The material 
parameters are shown in Table [3] The bulk modulus /fo is var- 
ied with respect to two values of the mean pressure /;>„, in the 
compressive strength relation (0.26 kPa and 1.3 kPa). During the 
impact a small region of the bottom of the dust sphere is com- 
pacted. Then the deformed sphere bounces off the target with 



reduced vel ocity Vf (see Fig. [TTIi . The latter effect was already 
observed in iGiittler et al.l (l2009bh and demonstrates the ability 
of our code and the implemented porosity model to simulate the 
elastic properties of the dust correctly. In this study we have dou- 
bled the resolution and determine the coefficient of restitution 
firest = ^7^o' depending on (see Fig.fTSTl. 

The bouncing calibration setup is equivalent with two cen- 
trally collidin g dust aggregates at 0.4 m s"', thus, our results are 
comparable to HeiBelmann et al.l (l2007h : Siest ~ 0.2 and « 95 % 
energy dissipation. 

Based on the results of the previous section, where = 
0.26 kPa turned out to be a good choice for the mean pressure, 
our results of this experiment favour a bulk modu lus A"o ^ 5 kPa. 
This value is close to - \ kPa computed by IWeidling et al.l 
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Fig. 15. Total cumulated volume over filling factor for drop ex- 
periments (crosses and triangles) and our 3D simulation with 
mean pressure - 0.26 kPa (solid line). The plot shows the cu- 
mulated mass fraction (normalised through the sphere volume) 
with a filling factor > (p over <p. Our simulation is in good agree- 
ment with the experimental findings. However, a higher amount 
of volume compressed to high filling factors leads to an almost 
constant deviation for (f> < 0.26. The slope is reproduced very 
well. 



Physical Quantity 



Symbol 



Value 



Unit 



Glass plate 



Bulk density'*' 


Po 


2540 


kgm 3 


Bulk modulus'*' 


Ko 


5 X 10' 


Pa 


Mumaghan exponent'*' 


n 


4 




Radius 


^ phite 


0.8 X 10-3 


m 


Thickness 


Opiate 


0.1 X 10-3 


m 


Dust sample 








Initial density 


p. 


300 


kgm"' 


Bulk density 


Ps 


2000 


kg m"3 


Reference density 


P'o 


300 


kg mr^ 


Initial filling factor 


0.15 




Bulk modulus 


Ko 


various 


Pa 


ODC mean pressure 


Pm 


260 and 1300 


Pa 


ODC max. filling factor 


<p2 


0.58 




ODC min. filling factor 


01 


0.12 




ODC slope 


A 


0.58 




Impact velocity 




0.2 


m s ' 


Radius 


r 


0.5 X 10-3 


m 



Table 3. Numerical parameters for the bouncing calibration 
setup. ODC stands for omni-directional compression relation 
('Eq.l211i. Qu antities m arked by (*) represent the parameters for 
sandstone in lMeloshl(fT989l) . which we adopt for glass here. 



(I2009h with a rough model. Our simulations yield a coefficient 
of restitution Siest - 0.19 (~ 96 % energy dissipation) for K^) - 
5.0 kPa, which is in excellent agreement with the experimental 
results. On the other hand for Ko = 500 kPa we find £,est - 0.09 
(«i 99 % energy dissipation), which is too far away from the 
reference. A high value for the bulk modulus Kq as it is given by 
the sound speed measurements is therefore excluded. 

Given the higher value p,,, - 1.3 kPa for the compres- 
sive strength curve Siest is raised for all choices of Kq. For 
Ko = 1.0 kPa it yields Siest ~ 0.7 and only 50 % of the 
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Fig. 16. Cross section through glass bead (red) and dust sample 
(light blue) at maximum intrusion for our 3D simulation (with 
Pm - 0.26 kPa, top) and the drop experiment (bottom). The 
colour indicates the spatially averaged filling factor The den- 
sity structure underneath the glass bead match very well. Even 
the slight compression along the tight intrusion channel can be 
reproduced. In the simulated plot a gap between glass bead and 
the most dense area is clearly visible. This is due to the smooth- 
ing of the sphere and was discussed in Sect. 14.21 The gap has 
roughly the size of one smoothing length h. 



energy is dissipated. On the other hand for Kq - 300 kPa we 
find firest ~ 0.13 equivalent to x 98% energy dissipation. For 
Pm = 1.3 kPa it is harder to dissipate energy by plastic deforma- 
tion. Hence, less energy is dissipated in the compaction process 
and the coefficient of restitution is generally higher 

This bouncing experiment fixes our choice for the bulk mod- 
ulus to Ko a: 5 kPa while maintaining - 0.26 kPa from the 
compaction experiment in the previous section. The next section 
will show that this choice is also consistent with the fragmenta- 
tion behaviour of the dust aggregates. 



5.3. Fragmentation 

Since the intended field of application of our calibrated SPH 
code will be the simulation of pre-planetesimal collisions, it is 
of major importance to calibrate and test the fragmentation be- 
haviour of the simulated material. For this reason we simulate 
the fragmentation experiment described already in the second 
part of Sect. [331 
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Fig. 18. Coefficient of restitution grest over bulk modulus Kq for 
values of the mean pressure /?,„ of the compressive strength re- 
lation. For the lower p,„ = 0.26 kPa (squares) more energy is 
dissipated by plastic deformation. For this reason the coefficient 
of restitution is significantly lowered compared to the higher 
Pm -1.3 kPa (crosses). Best agreement with the experimental 
reference of Siest = 0.2 (dashed line) is given for Ko x 5 kPa and 
Kq X! 20 kPa, respectively. 
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Fig. 19. Cumulative mass distribution of the fragments of a dust 
aggregate impacting on a glass plate for different values of Kq. 
For low fragment masses the shape of all simulated curves dif- 
fers from the experimental curve in Fig.|2]due to the limited res- 
olution of the experimental setup. An increase of Kq leads to an 
increase in the slope a of the power-law fit. The best agreement 
with the experimentally measured slope 0.67 was found for the 
simulation with Kq = 4.5 kPa. 



In our simulation a dust aggregate (0, — 0.35, 189296 par- 
ticles) hits a glass plate (188478 particles) from below with 
an impact velocity of vg = 8.4ms"'. The spatial resolution 
He = 8.0/im, h - 30.0/im) was chosen such that a single parti- 
cle has less than 5 x 10"* times the mass of the whole aggregate 
(6.8 X 10"^ kg) to resolve the same fragment masses as the ex- 
perimental reference. Gravitation was taken into account. A list 
of the other relevant parameters can be found in Table |4] 

We investigated the influence of the bulk modulus on the 
fragmentation behaviour. For this reason we varied Kq from 
3.0 kPa to 6.5 kPa. For comparison with the reference experi- 
ment (Fig. |2| we also plot our fragmentation data in a cumu- 
lative way (Fig. \T% and find also a good a agreement with the 
power-law distribution of Eq. |24] The simulation was evaluated 
after 0.8 ms simulated time. The mass of a fragment is given by 
the sum of the mass of the SPH particles belonging to it. We 
consider two fragments as separated when they are not linked by 
SPH particles that interact with each other, i.e. when the clos- 
est SPH particles of two fragments have a distance more than a 
sm oothing length . 

iGiittler et al.l (l2009bl) were facing the problem that using 
Kq - 300 kPa almost no fragmentation occurred (see Fig. 21 
in this reference). Their speculations about a modification of 
the shear strength in order to resolve this problem proved to 
be wrong. The quantity with the most impact on the fragmen- 
tation behaviour of the dust aggregate is its bulk modulus. In 
general it can be said that an increase of the bulk modulus leads 
to an increasing slope a of the fragment distribution, whereas 
the size of the largest fragment (normalised through the total 
mass of the projectile) fi roughly remains constant at ^ 20 % 
up to K() = 6.0 kPa. For higher Kq mainly small chunks and 
single SPH particles will burst off the aggregate, which is only 
compacted but does not fragment . This reproduces the situation 
described in lGuttler et"al] ( l2009bl) and Sect. [3] 

We calibrate for a, which is more sensitive to changes of Kq 
(see Table |5]). Given the measured value of a = 0.67 we find 
an excellent agreement with our simulation using Kq - 4.5 kPa, 



which yields a - 0.673 + 0.017. Moreover this simulation re- 
produces also the experimentally measured normalised mass of 
the largest fragment jj. - 0.22 to a very high accuracy (yi = 
0.234 + 0.007). The slight overestimation may be caused by the 
fact that the increase in filling factor is not taken into account in 
the analysis of the experimental data (see Sect.|3]l, whereas in the 
simulation it is. The fragment distributions for different Kq and 
the best fit for the power-law are shown in Fig. [19] Setup and out- 
come of the simulation are displayed in Fig.|20] As already indi- 
cated in Sect. l5.2l the choice of p,,, - 0.26 kPa and Kq = 4.5 kPa 
is consistent with the results from the compaction and bouncing 
experiments. The fragmentation experiment serves as a proof for 
the validity of these choices and also for the consistency of the 
underlying porosity model. 

In contrast to the experiments we find no material sticking 
to the glass plate. This is not possible due to the setup of the 
simulation. As in Sect. 14.41 we use artificial viscosity to separate 
glass and dust materials. This leads to an additional pressure on 
the dust material which prevents sticking. 

6. Discussion 

In this work we presented a smooth particle hydrodynamics 
(SPH) code equipped with extensions for conti nuum mechan- 
ics of solid bodies and an extended version of the ISironol (l2004l) 
porosity model. The code uses experimentally measured macro- 
scopic parameters such as tensile strength an d a static compres- 
sive strength relation. In lGuttler et al.l (l2009bl) this code was used 
to determine a relation for the shear strength and an estimate for 
the mean pressure p,„ for the dynamic compressive strength re- 
lation (Eq.|2TJ. The estimate was quite crude, though, due to the 
usage of 2D simulations only. 

This work profoundly investigated the numerical properties 
of the SPH code. Utilising the compaction calibration setup of 
Guttler et al. (2009b) as an example we determined an adequate 
size of the computational domain and adequate numerical and 
spatial resolutions. We have shown that the results for this setup 
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Fig. 20. Fragmentation sequence for f = ms (top left), t - 0.4 ms (top right), and t - 0.8 ms (bottom). The colour code indicates 
the filling factor A projectile consisting of dust particles (Sirono EOS, diameter 0.57 mm) hits a solid surface simulated by boundary 
glass particles (Murnaghan EOS, diameter 1.6 mm, thickness 0.04 mm) with an impact velocity of vq = 8.4 m i"'. A bulk modulus 
of Kq - 4.5 kPa was used for this simulation. While many small pieces, often single SPH particles, chip off, most of the aggregate is 
compacted to its maximum filling factor 0niax - 0.58 and fractures. The bottom layer of the impacting sphere is hardly compacted. 
This uncompacted layer is still visible on the fragments. (An animation of this figure is available in the online journal.) 



Physical Quantity 


Symbol 


Value 


Unit 


Glass plate 


Bulk density'*' 


Po 


2540 


kgm ^ 


Bulk modulus**' 


Ko 


5 X 10' 


Pa 


Murnaghan exponent'*' 


n 


4 




Radius 


f plate 


0.8 X 10-3 


m 


Thickness 


'Opiate 


0.04 X 10"' 


m 


Dust sample 








Initial density 


Pi 


300 


kg m" ' 


Bulk density 


Ps 


2000 


kgm"'' 


Reference density 


Po 


700 


kgm"^ 


Initial filling factor 




0.35 




Bulk modulus 


Ko 


various 


Pa 


ODC mean pressure 


Pm 


260 


Pa 


ODC max. filling factor 


(p2 


0.58 




ODC min. filling factor 


01 


0.12 




ODC slope 


A 


0.58 




Impact velocity 


vo 


8.4 


m s"' 


Radius 


r 


0.285 X 10-3 


m 



Table 4. Numerical parameters for the fragmentation calibra- 
tion setup. ODC stands for omni-directional compression rela- 
tion (Eq.l2TI). Q uantities marke d by (*) represent the parameters 
for sandstone in lMeloshl(ll989h . which we adopt for glass here. 



^0 [kPa] 


Slope 


a 


Norm, largest fragment fi 


3.00 


0.361 ± 


0.004 


0.200 


+ 


0.008 


3.50 


0.429 ± 


0.002 


0.172 


+ 


0.002 


4.00 


0.518 ± 


0.011 


0.230 


+ 


0.009 


4.25 


0.523 ± 


0.006 


0.194 


+ 


0.004 


4.50 


0.673 ± 


0.017 


0.234 


+ 


0.007 


4.75 


0.834 ± 


0.025 


0.196 


+ 


0.005 


5.00 


0.832 ± 


0.063 


0.198 


+ 


0.011 


5.50 


0.836 ± 


0.052 


0.220 


+ 


0.010 


6.00 


2.027 ± 


0.121 


0.171 


+ 


0.002 


6.50 


0.910 ± 


0.053 


0.390 


+ 


0.013 



Table 5. Results from the fragmentation calibration setup. The 
slope a of the power-law is increasing with increasing bulk mod- 
ulus Kq. Remarkably, the size of the normalised biggest fragment 
remains nearly constant around jj a; 0.2 for Kq < 6.0 kPa. 



are converging for higher spatial resolutions. Boundary condi- 
tions are necessary for all calibration setups presented in this 
work. Their treatment, which is a difficult issue in SPH, was re- 
solved by using boundary particles with vanishing acceleration 
at every time step. The dissipative properties of the artificial vis- 
cosity and its role for the stability of the simulation were investi- 
gated. Artificial viscosity was used to separate dust and glass ma- 
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terial which are highly different in the "stiffness" of their equa- 
tions of state. 

We investigated the crucial differences between 2D and 3D 
compaction calibration setup and proved the validity of the cor- 
rection factor for the intrusi on depth (Eg. [26ll tha t was already 
used without confirmation in lGiittler et al.l (l2009bl) . A major dif- 
ference between experiments and simulations, a compaction that 
went too far down in the dust sample and along with this a com- 
paction of too much volume, was resolved by using the 3D setup. 

Using a series of 2D simulations of the compaction calibra- 
tion setup we predicted a new, more accurate, value for the mean 
pressure p„, of the dynamic compressive strength relation. A 3D 
simulation with this value p,,, = 0.26 kPa was carried out. The 
results were compared with the experimental referen ce using the 
same benchmark features as in lGiittler et al.l (l2009bh . We found 
a goo d agreement. 

In iGiittler et al.l (l2009bl) it was already demonstrated that 
our code is in principle capable of simulating the bouncing of 
dust aggregates. With the bouncing calibration setup we now 
investigated the ability of the SPH code of quantitatively sim- 
ulating the elastic properties of the dust and the energy dissi- 
pation by compaction. We found that the bulk modulus as well 
as the compressive strength relation have significant influence. 
For smaller values of p„, (low compressive strength) and higher 
values of the bulk modulus more energy is dissipated. Using 
p,„ - 0.26 kPa from the compaction calibration setup we de- 
duced that the uncompressed dust samples {<p ^ 0.15) have a 
bulk modulus Kq x 5 kPa. With this result we were able to de- 
cide between two differi ng experimental values in favour of the 
indirect determination bv lWeidling et al.l (|2009|) . 

An important application of this code will be pre- 
planetesimal collisions. Therefore we also tested the code for 
its ability to simulate fragmentation of dust aggregates quantita- 
tively correct. For this we used a fragmentation calibration setup 
and identified the bulk modulus as the quantity with the most in- 
fluence on the fragment distribution. Again using p,„ = 0.26 kPa 
from the compaction calibration setup we found the best agree- 
ment with the empirical reference for the bulk modulus Kq - 
4.5 kPa. Remarkably this is consistent with the findings from the 
bouncing calibration setup, which represents a test for a totally 
different behaviour. The prob lem of almost no fragm entation for 
Ko = 300 kPa described in IGiittler et al.1 (l2009bl) was hereby 
traced back to a wrong assumption for the bulk modulus. 

7. Conclusion 

ISchafer et al.l (l2007l) have shown that the outcome of pre- 
planetesimal collisions crucially depend on their material prop- 
erties. If the issue of collisional growth is to be investigated by 
computer simulations, they suggest an implementation of experi- 
mentally measured material parameters and thorough calibration 
and comparison with laboratory experiments. 

The SPH code presented in this work complies with their 
requirements . Based on the experimental preparatory work by 
IGiittler et al.l (l2009bl) this code has been successfully tested with 
three kinds of calibration setups for the correct simulation of 
compaction, bouncing, and fragmentation. 

We conclude that we have developed a tool that features a 
sufficient accuracy for the investigation of the outcome of pre- 
planetesimal collisions in parameter ranges that are not accessi- 
ble to laboratory experiments. This is based on the assumption 
that the macroscopic properties calibrated for in this work do not 
change with increasing size. However, due to the thermodynam- 
ically simple porosity model used in the SPH code the area of 



application is restricted to a certain range of collision velocities. 
Collisions where shock propagation cannot be neglected are out- 
side the physically meaningful limit of this model. Changes of 
the mechanical properties with temperature and other thermody- 
namical effects like sintering and melting cannot be simulated. 
An extension of the porosity model and its implementation and 
calibration are necessary to broaden the parameter range (e.g. 
supersonic impacts) and to consider a realistic environment for 
pre-planetesimal collisions in the protoplanetary disc. 

Nevertheless, within its area of application our code will be 
used to produce a "catalogue of collisions". The influence of ob- 
ject sizes, porosity, relative velocity, rotation, and impact param- 
eter on the fragment distribution of pre-planetesimal collisions 
will be investigated in future work. 
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